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PsJ , Abstract 

■ Model applicability to core-scale solute transport is evaluated using breakthrough 

j^ . data from column experiments conducted with conservative tracers tritium (^H) and 

sodium-22 (^^Na), and the retarding solute uranium-232 (^^^U). The three models 
_. considered are single-porosity, double-porosity with single-rate mobile-immobile mass- 

i-^ . exchange, and the multirate model, which is a deterministic model that admits the 

c/3 I statistics of a random mobile-immobile mass-exchange rate coefficient. The experi- 

^^ ' ments were conducted on intact Culebra Dolomite core samples. Previously, data were 

£^' analyzed using single- and double-porosity models although the Culebra Dolomite is 

t— I ! known to possess multiple types and scales of porosity, and to exhibit multirate mobile- 

Qh| immobile-domain mass transfer characteristics at field scale. The data are reanalyzed 

here and null-space Monte Carlo analysis is used to facilitate objective model selec- 
tion. Prediction (or residual) bias is adopted as a measure of the model structural 
error. The analysis clearly shows single- and double-porosity models are structurally 
^T) • deficient, yielding late-time residual bias that grows with time. On the other hand, 

C^^ . the multirate model yields unbiased predictions consistent with the late-time —5/2 

slope diagnostic of multirate mass transfer. The analysis indicates the multirate model 
is better suited to describing core-scale solute breakthrough in the Culebra Dolomite 
than the other two models. 
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^ . 1 Introduction 

_^_' During the last 30 years, si Knificant effort has be e n expended to understan d contaminant 



transport in fractured rock [Huyakorn et al.l . Il983l . ISun and Buscheckl . l2003l | due in part to 



the necessity to evaluate site suitability for nuclear waste disposal. Contaminant transport 
in fractured rock is of common concern to regulators and stakeholders at nuclear waste 
disposal sites because off-site contaminant migration could impact groundwater resources. 
Modeling contaminant transport in fractured rock is challenging due to the complex and 
inherently heterogeneous nature of the transport domain, and the multitude of physical and 
chemical processes controlling contaminant interaction with the host rock. This has led to 



a deve lopment of several potentially competing conceptualizat ions of the transport environ 



van Genuchten and Wagenetl . Il989l . IZheng et all |2010|. Model selection is typically 



ment 

based on subjective expert judgment. Hence, there is a need for objective criteria for se- 
lecting physically-based models that best describe observed transport behavior and provide 
minimal predictive uncertainty. 

In this work, we present a criterion for selecting between competing models for describ- 
ing transport at the core scale. Three models are considered: the single-porosity model; the 
traditional double-porosity model with single-rate mob i le-im mobile domain mass exchange 



van Genuchten and Wagenetl . ll989l . lGamerdinger et al.l . ll990j . and; a double-porosity model 



with multiple rates of mobile-immobile-domain mass ex change controlled by a random mass 



transfer coefficient [Haggerty and Gorelickl . 1 19951 . Il998| . We refer to the traditional double 



porosity model as simply the double-porosity model, and to the model with multiple rates of 



mass exch ange as the rnultirat e model following iHaggerty and Gorelickl [l995J . lHaggerty et al. 



(2000l | and iMeigs et al.l 2000| . In the multirate model, the mass transfer coefficient is a ran- 
dom variable, not a single deterministic parameter. This conceptualization reflects spatial, 
not temporal, variability (due to heterogeneity, i.e., multiple types and scales of porosity). 
The probability density function of the transfer coefficient gives the probability that a mobile- 
immobile interface (assumed to be randomly distributed in space) , encountered by a particle 
along its trajectory through the transport domain, has a particular mass transfer coefficient 
value. 

The three mo dels are used to ana lyze breakthrough data collected in core-scale labora- 



tory experiments Lucero et al.l . Il998| using conservative tracers tritium (^H) and sodium-22 
(^^Na), and the retarding tracer uranium-232 (^^^U). The experiments analyzed herein were 
performed on a rock core collected from a formation known to exhibit multiple types and 
scales of rock-matrix poro sity. Previous analy sis of the experimental data with single- and 
double-porosity models by iLucero et al.l |l998| yielded poor model fits to these data due to 
the inability of the two models to describe the long-tailing behavior of conservative solutes. 
The multirate model has been shown t o properly describe this behavior in breakthrough 



data obtained in fleld -scale tracer tests [Meigs and Beauheiml . |2001| . IHaggerty et al.l . 12001 



McKenna et al.l . l200l| . It is applied herein for the flrst time to core-scale breakthrough data 
to demonstrate multirate mass-transfer effects are observable at this scale. 

Null-space Monte Carlo analysis (NSMC) is used to evaluate model prediction uncertainty 
for each of the three mod el based on br e akthrough data. It yields rnultiple sets of parameters 



that calibrate the model Tonkin et al.l. l2007l. iTonkin and Dohertyl. l2009l. iJames et al.l. 12009 



Gallagher and Dohertyl . |2007|, leading to multiple realizations of model flts to data at pa- 
rameter estimation optimality. By prediction uncertainty we mean the variance and bias of 
the ensemble of these model-prediction realizations relative to observed behavior. Variance 
describes the scatter of realizations about mean behavior, while the residuals bias associated 
with each data point at optimality over all NSMC realizations provides a measure of the 
systematic departure of predicted from observed behavior. This work presents the flrst use 
of residual bias in the solute transport literature as a criterion for model selection. 



2 The multirate transport model 

The multirate model is based on the traditional double porosity model where the transport 
domain is conceptualized as comprising two overlapping continua, namely the mobile (ad- 
vective or fracture porosity) and immobile (diffusion-dominated matrix porosity) domains. 
Unlike the traditional double porosity model where a single deterministic constant is used 
to characterize mobile-immobile-domain mass exchange, a random variable is used in the 
multirate model. Using this conceptual app roach, the governing equation for t ransport of 



a sorbing radionuclide in the mobile domain |Haggerty and Gorelickl . 119951 . Il998l | is given in 
nondimensional form 

where C = c/Cc, Cjm = Qm/C'c, X = x/Lc, T = t/Tc, c and cim are mobile- and immobile- 
phase solute concentrations [M L~^], x and t are space-time coordinates, Cc, Lc, and Tc 
are characteristic concentration, length, and time. Ad = XTc, A is the first-order radioactive 
decay constant [T^^], ud = ojTc is the dimensionless first-order mass-transfer rate coefficient 
(Damkohler-I number), [3{ujd) = (3tp{^d) is the rock matrix point capacity ratio, (3t = 
4>imRim/(pmRm is the dimeusiouless rock-matrix total capacity ratio, p{ojd) is the probability 
density function (pdf) of ojd, Pe = LcJcxl is the Peclet number, a^ [L] is the longitudinal 
dispersivity, 0m and ^im are the mobile- and immobile-domain porosities, and i?m and R\^ 
are the mobile- and immobile-domain retardation factors. 

The dimensionless governing equation for immobile domain transport is 

— — — = UJd{C — Cim) — A/jCim, (2) 

ol 

the lumped-parameter formulation of immobile-domain mass transport. 
The transport equations are solved subject to the initial condition 

c(x,r = o) = am(r = o) = Co, (3) 

indicating initial equilibrium between mobile and immobile-domain concentrations. The 
boundary condition at X = is 



UdX 



X=0 



an,{T), (4) 



where Cinj is a normalized injection concentration and A [L] and B are parameters to specify 
the X = boundary condition type {A = and B = 1 correspond to a Dirichlet boundary 
condition, while A = —D/v and B = 1 correspond to a Robin boundary condition). The 
downstream boundary condition is 

lim f-L^ + c] =0, (5) 

x^oo V Pe 9X / ' ^ ' 



indicating zero solute flux infinitely far downstream. 

The solution to dl])-® is obtained on a semi-infinite don iain < X < oo as a simpli- 



ficati on and limiting case of the finite domain considered by iHaggerty and Gorelickl |1995 
It is given by 



1998 



CiX) 



^inj ^ -Dtvp 



^uX 



a 



v 



(6) 



B + uA/L^ 

where the overbar indicates the Laplace transform, s is the Laplace transform parameter, u = 
(l - Vl + 4/i/P) P/2, Cp = Co/{s + Xn), MXd) = {s + Xd)MXd), MXd) = 1 + /3t^(Ab), 
and 

mo) = r :f '7' d... (7) 

Jo s + Xd + ujd 



The function g{X£)) is the Laplace transformed memory function of IHaggerty et al.l 20001 . 
For single porosity g{Xr)) = 0, whereas for double porosity with single-rate mass transfer 
QJXn) = 0->n/{s + An + Wd)- The inverse Laplace transform of (jS]) is obtained using the 



de Hoog et al.l |1982| algorithm. For all results reported herein, Cc = Cinj is the injection 



concentration, Lc is core length, and Tc = Lc/vr, where vr = v/Rm and v is the average 
linear velocity [L T~^]. 



2.1 Mass-Transfer Coefficient Distribution 

To evaluate the memory kernel (?(Ad) numerically, the probability density function p{ud) 
must be specified. All valid probability density functions are admissible in the computation 
of the memory function, incl uding single-paramete r distributions such as the power-law used 



by IHaggerty et al.l |2000| and lSchumer et al.l [2003|. However, single-parameter distributions 



may not lead to improved multirate model predictions of breakthrough behavior compared to 
the single-rate mass transfer model. Here, without loss of generality, we use the lognormal 
distribution because several key geologica l properties appear to approxima tely follow this 



distribution iHaggertv and Gorelickl . |1998| . inclu ding hydraulic condu ctivity [Neumanl . 11982 



Hoeksema and Kitanidisl . Il985| and grain size Buchan et al.l . Il993| . Other equally valid 



examples of distributions that have been used in the lit erature to characteriz e the mobile 



Haggerty et al.l [2000|. Using any 



immobile mass transfer coefficients are summarized in 

of these models with two or more parameters, would likely yield multirate models that 

outperform the single-porosity and single-rate double-porosity models. 

Th e stan dard two-parameter lognormal distribution for ud G [0, oo) was used by IHaggerty and Gorelick 



1995Ul998| . For the case where physical bounds exist ud G [wD.min, i^D.max], it may be more 
appropriate to use the random variable ud = (I/^d — l/wD.max)" — i^D.min, where uJD,mm 
and ujD^raax are the minimum and maximum physically allowable ud values. The pdf of cud 
is 



P(^i?) 



1 



UDd'\'2n 



exp 



log(^D) - A 
aV2 



(8) 



where fi and a are the mean and standard deviation oilog{ujo)- This is the lower- and upper- 
tail-truncated lognormal distribution, which is a more plausible distribution when there are 



physical limits on permissible Damkohler-I numbers. These physical limits may be estimated 
from data. In the limit as a;_D,mm — > and w^^max -^ co, '^d -^ <^d, and p{ojd) degenerates 
to the standar d two-parameter lognor r nal di stribution. We set ud G [0, 1000], loosely based 
on the work of iHaggerty and Gorelickl 1995| . where ud = 100 was suggested as the limit for 
significant multirate mass transfer. 



3 Application to Core-scale Breakthrough Data 

The data considered here were collected in a series of column experiments con ducted on five 



intact cores (denoted A through E) of the Culebra Dolomite as reported by iLucero et al. 



1998 . The Culebra Dolomite member of the Rustler formation of the Permian Basin in 



southeastern New Mexico is known to exhibit several categories and scales of porosity [Holt 



19971 ] including inter-crystalline, inter-particle, fracture, and vuggy porosities (Figured]). The 
multiple types and scales of porosity are also clearly observable in Culebra Dolomite cores 
(Figure [2]). The only breakthrough data analyzed in this work were collected on the B core 
for the conservative tracers '^H and ^^Na, and the retarding tracer ^'^^U. Core B, pictured 
in Figure [H was selected because its length-to-diameter ratio (50.9 cm to 14.5 cm) ratio 
was such that boundary effects can be neglected, thus permitting the use of the analytical 
solution developed for a ID semi-infinite (0 < x < oo) transport domain. Dry bulk density 
Phuik = 2 400 kg/m^ and tota l porosity 0t = 0.14 were determined by standard laboratory 



methods [Lucero et al.l . 119981 ]. Additional d etails on experiment setup, solute injection, flow 
rates, and effluent analysis, are available in lLucero et al.l 19981 ] and are not repeated here. 

Figure [3] shows normalized concentrations plotted against pore volume (PV) computed 
using 0T- Solute injection pulses were longer in duration for tests shown in Figure [3](b) than 
for those in Figure [3](a). Plotting data on a log- log scale as in Figure [3](b) clearly shows that 
the effluent was not collected for a sufflciently long time to completely reveal the late-time 
tracer behavior. A long breakthrough tail is characteristic of mobile-immobile-domain mass 
transfer for conservative tracers. Despite this shortcoming, the data can be used to assess 
the performance of the three models. The data in Figure [3t^a) show early breakthrough 



for both conservative tracers [Lucero et al.l . 119981 ]. suggesting the occurrence of preferential 



flow in an advective porosity that is significantly smaller than the total core porosity 0t- 
Breakthrough data for 232-[j g^j,g shown in Figure|l](^^Na data from the same test are included 
for comparison). ^32^ breakthrough clearly occurs much later than ^^Na because the former 
sorbs onto the Culebra Dolomite. Peak ^32^ concentration arrival occurs around 1 PV, about 



four times later than ^^Na. Using the single-porosity model, iLucero et al.l 1998] estimated 
the 232-[j retardation factor to be 4.5 and 3.7, from B3 and B7 data, respectively. For the 
dual-porosity model, they obtained mobile- and immobile-zone retardation factor values of 
{/2in = 1.14, _Rim = 65.4} and {R^ = 4.35, -Rim = 1.00}, from B3 and B7 data, respectively. 
The value of -Rim = 65.4 appears to be an error in recording the estimated value. 



3.1 Parameter Estimation 

To estimate model parameters we let Cobs be the breakthrough data vector, Ccai(^) the 
model-calculated concentrations vector, and the vector of estimated model parameters. 
For ^H and ^^Na, 6 = (^m, aL,/i, a, tinj), whereas for ^^^U, 6 = (_Rin, _Rim,/i, a). Injection 
pulse concentration (cinj) was fixed for tests Bl, B2, B3, and B7, but was estimated for 
tests B4, B5, and B8. Increased test durations for B4, B5, and B8 made it more difficult to 
maintain constant injection concentrations over pr olonged test periods , resulting in injection 



concentrations that varied appreciably with time [Lucero et al.l . |1998| . Since this temporal 
variability is not incorporated explicitly into the solution, and its functional form in unknown, 
the injection concentrations for tests B4, B5, and B8 are treated as unknown constants and 
are estimated from breakthrough data. Initial concentration (cq) was fixed for all tests and 
was determined from effiuent concentration values measured prior to solute injection. The 
truncated lognormal distribution {ud G [0, 1000]) was used to describe the mass-transfer 
coefficient distribution. The advective porosity {4>m), dispersivity, and the injected pulse 
(tinj) duration were estimated with the multirate model for ^^Na data and used as fixed input 
parameters when estimating the retardation factor and ud distribution parameters from 232"u 
data. Distribution parameters were also estimated for ^sa-jj because ud is a function of the 
tracer-specific molecular diffusion coefficient. 

We examine model sensitivity coefficients to determine whether all model parameters are 
estimable from available data. Sensitivity coefficients are derivatives of model-predicted ef- 
fluent concentrations with respect to model parameters, which are elements of the Jacobian 
matrix (J). They provide a measure of parameter identiflability, beca use the determinant 



of_JjJ must be sufficiently larger than zero to be estimable from data Ozisik and Orlande 
2000|. Small sensitivity coefficients imply |J^J| ~ and the inverse problem is ill condi- 



tioned. Here, sensitivity coefficients were estimated with PEST using central differences, 
and their variation with time is shown in Figure |5]for (a) short (B2) and (b) long (B4) solute 
injection pulses. The sensitivities are sufficiently larger than zero to permit estimation of all 
parameters from breakthrough data. The coefficients are also linearly independent for much 
of the time data were collected. Apparent linear dependence is restricted to late-time data, 
implying parameters cannot be uniquely estimated solely from late-time data. The param- 
eter sensitivity curves obtained in both short- and long-pulse injection tests show a weak 
symmetry between two opposite-sign branches associated with arrival and elution tracer 
breakthrough waves. Absolute values of sensitivity coefficients are largest when measured 
concentrations are changing most rapidly. Variation of sensitivity coefficients with time for 
retarding tracer 232-[j jj^ ^gg^ g3 g^j,g shown in Figure |6l These are also sufficiently larger than 
zero indicating that parameters, including i?m and Rim, are estimable from breakthrough 
data. 



Parameter estimation was performed using PEST |Dohertyl . |2010| . The optimal vector 



of model parameters (^opt) was obtained by minimizing the sum of squared residuals, 

$(0) = eieyeiO), (9) 

where e = Cobs — Ccai(^) is the vector of residuals. PEST uses the Levenberg-Marquardt non- 



linear optimization algorithm Marquardtl . Il963| . Parameter estimates and multirate model 



fits to data are compared to those obtained using single- and double-porosity models. Param- 
eter values obtained by inverting ^H and ^^Na breakthrough data with the all three models 
are summarized in Table [U parameters estimated from ^■^^ U data are in Table [2j Because tinj 



was not reported in the original study [Lucero et al.l . Il998| , it was estimated from data. The 



ujo column also includes the mean ((wd)) and variance (a^ ) of the Damkohler-I number 
determined from the estimated values of fi and a. The last row of Table [T] shows estimated 
model parameters from simultaneous inversion of B4, B5, and B8 tracer-test breakthrough 
data. Parameter estimates are comparable to those from individual tests, even though the 
three tests were conducted with flow rates ranging over an order of magnitude (0.05, 0.1, 
and 0.5 ml/min). This indicates minimal model structural error with regard to simulating 
average pore-water velocity. 

Model fits to data for parameter values listed in Table [T] are shown in Figure [7] (BI- 
BS, B7) and Figure [8] (B4, B5, B8) for ^H and ^^Na. Figures are in pairs of (a) linear or 
semi-log (concentration on linear scale) and (b) log-log plots, to illustrate how models match 
data over multiple time scales and over several concentration orders of magnitude. The two 
plotting scales are complementary because an apparently good model fit on a semi-log or 
linear plot may reveal a po or fit on log- log sca le, and vice versa. Model-fit results for 232]j 



data are shown in Figure [H iLucero et al.l 1998| parameter estimates are comparable to those 
obtained here using single- and double-porosity models, but they did not estimate tmj. 

Parameter estimation using the multirate model yielded improved model fits to break- 
through data compared to those obtained using single- and double-porosity models (see B^ 
values in Tabled]). Mobile-domain porosity values (0m) estimated with single- and double- 
porosity models were comparable (means of 0.069 and 0.065, respectively), but were appre- 
ciably larger than those obtained using the multirate model (mean of 0.045). Dispersivity 
(ol) values were consistently largest for the single-porosity model (mean of 12.1 cm) and 
smallest for the multirate model (mean of 3.76 cm) for all tests. Table [1] shows there is signif- 
icantly more variability in a^ estimated using the single-porosity model than those obtained 
using the double-porosity and multirate models (standard deviations of 4.2 cm, 2.4 cm, and 
2.3 cm, respectively). The Damkohler-I numbers estimated with the double-porosity model 
appear closer (though not equal) to the geometric mean {{ujD)g = e^) of the multirate model 
than to the mean {{ujd) = e'^"'"'^ /^). Results show absolute values of /i and a for the ^H 
tracer test (Bl) are smaller than those obtained with the tracer ^^Na. With exception of 
B7, the ^^Na tests yielded consistent values of /i and a with |/i| > 1.0 and a ~ 1.9. Those 
obtained from ^^^-y ^jg^^g^ (Table [2]) are significantly different. 

For the non-conservative tracer ^^^U, 0m and a^ were estimated with the ^^Na tracer from 
the same experiment, because these parameters are intrinsic transport medium properties. 
Estimated retardation factors from tests B3 and B7 are listed in Table [2J For test B3, fitting 
the multirate and double-porosity models to data yields -Rm values appreciably smaller than 
the value obtained with single-porosity model. This is because retardation is distributed 
between the mobile and immobile domains in the former two models. It is surprising to find 
the multirate model i?im in test B3 is significantly larger than the double-porosity model 



i?im. Intuitively, one would expect results similar to those obtained from test B7, because 
delayed breakthrough is partly due to matrix mass transfer and partly due to solid-phase 
sorption. In addition, the retardation factors, -Rm and Rim, estimated with the double- 
porosity and multirate models showed significant differences between test B3 and B7. These 
two results may be attributable to interplay between multirate mass-transfer and nonlinear 
sorption kinetics, where retardation is concentration dependent. The models all assume 
linear instantaneous sorption, variability in retardation factors between tests B3 and B7 
may be an artifact of inherent model deficiency to account for nonlinear sorption kinetics. 
^^^U column tests tests B3, B6 (not discussed here) and B7 were performed serially on the 
same core. B3 had the lowest initial relative ^32^ concentration with Co/cjnj ~ 2 x 10~^, 
while for B7 Co/cinj — 10^^. B7 was performed after the core had already been conditioned 
with 232ij fj-om the previous two tests. These initial concentration differences are expected 
to affect the estimated retardation factors in the presence of nonlinear sorption kinetics. 

3.2 Predictive Analysis 

All models approximate a complex reality, and the discrepancy between reality and mathe- 
matical models is commonly referred to as model structural error. It is a measure of model 
deficiencies that lead to prediction errors even when the models are supplied with optimal 
input param eters. Structural error cann ot be attributed to measurement errors inherent in 
observations Doherty and Welterl . l201Cll | and typically decreases as models become more real- 



istic with increased understanding of underlying causal mechanisms of processes. A measure 
of structural error would thus provide an objective criterion for model selection. 

Predictive uncertainty analysis presented here is used to demonstrate the structural defi- 
ciency of the single- and double-porosity models, and how this deficiency leads to increased 
model prediction error. The analysis was undertaken with PEST for test B 8. Details for 
conducting a PEST predictive uncertainty an a lysis can be found elsewhere 
20091 . iTonkin and Dohertvl . l2009l . iTonkin et all . l2007l . ICallagher and Dohertv . 



James et al. 



20071. Using 



parameter values at optimality (Table [T]) and the associated covariance matrix, 500 random 
parameter sets were generated and projected onto the Jacobian matrix null space. No clear 
null space was found from the singular value decomposition of the Jacobian matrix, therefore 
we assumed the null space to be a single dimension in these low-dimensional (< 6) models. 
Model predictions computed beyond the last observation based on the 500 parameter sets 
generated in this manner are shown in the left column of Figure (TU] for (a) single-porosity, 
(c) double-porosity, and (e) multirate models. They show significant model prediction un- 
certainty for the single-porosity model, and only moderate uncertainty for the other two 
models. Using these parameter sets projected onto the null space as initial guesses, further 
minimization of $ was undertaken, using the Jacobian matrix associated with the calibrated 
state. Using the value of $ at optimality (^opt), the 500 null- space-projected parameter sets 
were processed with PEST to minimize the objective function such that $ < 2$opt- Pre- 
dictions associated with the re-calibrated parameter sets are shown in the second column of 
Figure [To] for (b) single-porosity, (d) double-porosity, and (f) multirate models. As would be 
expected, post re-calibration model predictions for all three models show a marked decrease 



in model prediction uncertainty from the pre-calibration predictions. The late-time —3/2 
and —5 /2 slope lines a r e incl uded, which are diagnostic of double-porosity and multirate 



models jHaggerty et al.l . |2000| |. Clearly, the model behavior projected beyond the time of 
the last observation follows the —3/2 slope for the dual-porosity model, and the —5/2 slope 
for the multirate mass transfer model. 

Re-calibration single-porosity model projections show significant underestimation of late- 
time observations. Dual-porosity model predictions are skewed toward overestimating the 
late-time observations. Multirate model projections are uniformly centered about the data 
and are consistent with the observed trend of the elution curve. Figure [11] shows histograms 
of residuals associated with the three models plotted at (a) t = 4.1 and (b) t = 4.7 days. 
Whereas the residuals computed at t = 4.7 days with the multirate model have zero bias, 
those of the double- and single-porosity models show clear bias to negative (concentration 
overestimation) and positive (underestimation) values. Only the multirate model shows 
minimal bias about the observed late-time data, even though its ensemble of predictions 
has comparable spread (variance) to those of the double-porosity model beyond the last 
observation. The residual bias signifies model structural error associated with single- and 
dual-porosity models. Comparing results in Figure [TT] (a) and (b) shows residual bias and 
single- and double-porosity model structural error increase with time, while bias for the 
multirate model does not show appreciable change. At time t = 4.1 days, the dual-porosity 
model residuals have zero mean and are nearly coincident with the multirate model. However, 
at t = 4.7 days there is a growth in double-porosity model prediction bias. Prediction error 
due to model structural error increases with time. 

Figure [12] shows histograms of 500 calibrated multirate model parameter sets obtained 
from the posterior null-space Monte Carlo analysis described above. These dis tributions pro 



vide a measure of parameter estimation uncertainty. However, as indicated by [Keating et al. 



2010|, parameter sets obtained using null-space Monte Carlo analysis do not necessarily 
constitute a sample of the posterior density function of the parameters in the strict Bayesian 
sense. This is especially true with low-dimensional models (at most 6 parameters for the 
present case) for which a proper null space may not exist. This can be seen by comparing 
the posterior distribution obtained with the null-space Monte Carlo analysis with those ob- 
tained to a formal Ba yesian approach using the DiffeRential Evolution Adaptive Metropolis 



(DREAM) algorithm Vrugt et al.l . [20081 [2009a[ Jb| . For the problem considered here with 6 



parameters to be estimated from log-transformed concentrations, DREAM ran 6 different 
Markov chains, and after a burn-in period of about 35,000 model runs per chain, we obtained 
the parameter posterior distributions shown in Figure [13] DREAM required 300,000 total 
model runs. Clearly, the computational demands of formal Bayesian analysis with DREAM 



can be prohibitively high |Keating et al.l . [2010| . The parameter posterior distributions shown 



in Figure [T^ show the final 10,000 model runs. Normal distributions are included in the figure 
for comparison. The results show that posterior distributions obtained with DREAM have 
smaller variances and are more Gaussian than those obtained with the PEST posterior null- 
space Monte Carlo analysis. Whilst PEST results indicate greater variability in estimated 
parameter values that calibrate the model, DREAM results indicate that parameter estima- 



tion uncertainty is actually smaller. The low-dimensionality of the parameter space leads to 
an overestimation of parameter estimation uncertainty using null-space Monte-Carlo analy- 
sis. Thus, PEST-based parameter estimation uncertainty, obtained with null-space Monte 
Carlo analysis for a significantly lower computational cost, may be viewed as the upper 
bound of the true uncertainty computed with DREAM, for cases like the low-dimensional 
models used here. 



3.3 Statistical Model Selection 

For a given number of obervations, as models become more realistic, the increase in model 
complexity and the number of parameters leads to increased parameter estimation uncer- 
tainty because the number of observations available per estimated parameter decreases. In 
the present case, model complexity and the number of parameters increase from the single- 
porosity to the multirate model, but the respective model parameters are estimated with the 
same number of observations. Hence, statistical criteria that account for decreased informa- 
tion content due to increased model complexity may be used to augment model selection 
based on structural error evaluation. The corrected Akaike Information Criterion, AIC,. 



Hurvich and Tsai. 1997. Anderson and Burnham. 1999, Poeter and Anderson. 20051 is used 



here for this purpose 



AIC, = 2n 



log((Te) + 



n 



(10) 



where n is the number of observations, k is the number of estimated parameters, and Ue is 
the standard deviation of residuals at optimality. The first term typically decreases as model 
complexity increases, representing improved model fit to data, while the second penalty term 
increases. Because AIC^ is a relative measure, it is preferable to use differentials of AICc 
Posada and Buckleyl . 120041], denoted AAICc, over all the three models under consideration. 
For the i^^ model, AAICc,i = AICc,j — minAICc, where min AICc is the smallest AICc value 
among all models for this dataset. The AICc are computed using PEST and AAICc are 
listed in Table [H The minimum AICc corresponds to the multirate model, except in test 
B5, where the it corresponds to the double-porosity model. Clearly the relative AICc values 
confirm the results of predictive analysis that the multirate model is better suited than the 
other two models to describing transport in the Culebra Dolomite core. 

For time series data with high autocorrelation, the penalty for model complexity is van- 
ishingly small when n ^> k and the AICc reduces to a ranking of the models by residual 
variance. This is only a problem however, when the increased number of observations does 
not singinficantly increase the information content of the observation about the estimated 
parameters. Hence, a separate optimization with PEST using only 30 of the original 269 
data in test B8 to determine whether the ranking of the three models with the AICc would 
change appreciably. The resulting model fits are shown in Figure [HI Basically the same 
results were obtained, with the multirate model outperforming the other two models. This 
is because the estimation variance is always smallest for the multirate model, and artificially 
reducing n only has a modest effect on the final outcome. It should also be noted that a large 
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n allows one to better capture the variability in the data due to random measurement error, 
which are assumed to be Gaussian in minimization of the sum of squared residuals. Further, 
the number of parameters to be estimated increased only by 2 from the single-porosity to the 
multirate model, whereas the estimation variance changes by a factor of about 2 (7.6 x 10~^ 
to 3.2 X 10-^). 

The temporal structure of the residuals was examined to determine whether they show 
strong temporal autocorrelation. They are plotted in Figure [T51 It can be seen in the Figure 
that moderate autocorrelation is limited to very early-time. Additionally, in this early-time 
period, it can be seen that only the single-porosity model residuals show appreciable tem- 
poral autocorrelation, which decreases as one moves to the multirate model. The computed 
responses of the single-porosity model show strong departure from observed behavior. As 
can be seen in Figure [13 the residuals obtained with the multirate model for the long tests 
(B4, B5 & B8) show only moderate temporal autocorrelation (at early-time) and are mostly 
randomly distributed about zero. It should also be noted that the statistical rigor of DREAM 
does not depend on the distribution of the residuals but on the sampling of the parameter 
space for parameters that minimize the sum of squared residuals. 

4 Discussion and Conclusions 



We reanalyzed core -scale '^H and ^^Na breakthrough data from experiments conducted by 
Lucero et al.l |1998| on a C ulebra Dolomite core u s ing the sin gle-porosity, double-porosity, 
and the multirate model of Haggerty and Gorelickl . 119951 . Il998| on a semi-infinite domain to 
determine which of the models best describ es the observed breakthrough behavior. Previous 
analysis of these data by lLucero et al.l |1998| had suggested that the single-porosity model was 
sufficient to describe core-scale Culebra transport, a finding that was at odds with findings 



based on field-s c ale te sts conducted in the Culebra Dolomite formation [Meigs et al.l . [2000 



McKenna et al.l . |2001| . In the results presented herein, the multirate model yielded better 
model fits to the data and parameter values that differed significantly from those obtained 
with the single- and double-porosity models. The mobile-domain porosity and dispersivity 
values obtained with the multirate model were consistently lower than those obtained with 
the other two models because solute dispersion in the core is also accounted for by porosity 
variability encapsulated in the distribution parameters of the mobile/immobile domain mass- 
transfer coefficient. The smaller dispersivity obtained with the multirate model is more 
plausible than those obtained with the other models, considering the length scale of the 
experiments. 

Model-prediction uncertainty was evaluated using breakthrough data from test B8 and 
post-calibration null-space Monte Carlo analysis as implemented with PEST. The predic- 
tion uncertainty analysis revealed the presence of model structural error in the single- and 
double-porosity models as demonstrated by significant bias in the residuals of model predic- 
tions made with these models with optimal parameter values. The residual bias increased 
with time over the span of the elution curve where breakthrough data are available, show- 
ing increased departure of model predictions from the observed trend (—5/2 slope line) of 
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breakthrough data. The parameters associated with the null-space Monte Carlo predictive 
analysis may be viewed as samples from the posterior parameter distributions and were used 
to evaluate parameter estimation uncertainty. The posterior distributions estimated using 
null-space Monte Carlo analysis were compared to those obtained with the more rigorous 
Bayesian analysis in the DREAM algorithm. The comparison suggests that measures of 
parameter estimation uncertainty obtained with null-space Monte Carlo may be treated as 
upper bounds of the true posterior distributions, particularly for low dimensional models 
where a true null space may not exist. 

The analysis presented herein clearly shows the residual bias associated with the single- 
and double-porosity models increases with time indicating increasing systematic departure of 
predicted from observed behavior due to the inherent structural deficiencies of these models. 
The multirate model residuals, however, maintain minimal bias with time, indicating low 
model structural error. Although the predictions with the double-porosity and multirate 
models beyond the last observation have comparable variance, only the residuals of the 
multirate model have zero bias. These results show that the multirate model is the most 
appropriate of the three models for describing solute breakthrough behavior in Culebra 
core even though the three models yield parameters with comparable variances of posterior 
distributions. This finding was further confirmed using statistical model selection using the 
differential AICc where the AICc value was typically smallest for the multirate model. The 
one test where the double-porosity model yielded the smallest differential AICc value, the 
value associated with the multirate model was only marginally larger (0.5%). More elution 
data would be needed to resolve this minor departure from the norm given that the two 
models predict disparate long-term tailing behaviors. 
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Figure 1: Different types and scales of Culebra Dolomite porosity. Anhydrite and mudstone 
of adjacent Rustler members act as confining layers. 
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Figure 2: Culebra Dolomite horizontal core B showing vuggy porosity, fractures, and vug- 
filling calcite. Foreground grid marks are inches. 
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Figure 3: Normalized concentrations plotted against PV for (a) short-injection-pulse tests 
Bl, B2, B3, and B7, and (b) long constant-concentration-injection tests B4, B5, and B8. 
Vertical line marks one PV calculated using total porosity. Q^ in (b) is volume fiow rate in 
ml/min. 
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Figure 4: Concentrations plotted against PV for ^^Na (conservative) and ^ss-q (retarding) in 
tests B3 (a) and B7 (b). Vertical line marks one PV calculated using total porosity. 
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Figure 5: Breakthrough concentration sensitivities to estimated multirate model parameters 
for (a) short- (B2) and (b) long-pulse (B4) ^^Na tests. Concentration data are included for 
reference (those in (b) are scaled by 0.5). 
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Figure 6: Breakthrough concentration sensitivities to estimated multirate model parameters 
for 232U in test B3. 
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Figure 7: Model fits to core B ^^Na breakthrough data for short-pulse tests (Bl, B2, B3, 
and B7) on (a) semi-log and (b) log-log scales. 
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Figure 8: Model fits to core B ^^Na breakthrough data for long-pulse tests (B4, B5, and B^ 
on (a) linear and (b) log-log scales. 
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Figure 9: Model fits to ^sz-y breakthrough data from tests B3 and B7 on semi-log plots. 
Parameters shown in plots are for the multirate model. 
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Figure 10: Model prediction uncertainty evaluated using posterior Monte Carlo analysis on 
B8 data with the (a,b) single-porosity, (c,d) double-porosity and (e,f) the multirate models 
(left and right columns represent before and after re-calibration, respectively). Double poros- 
ity model predictions approach —3/2 slope while the multirate model predictions approach 
—5/2 slope. 
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Figure 11: Residual histograms computed at (a) t = 4.1 days and (b) t = 4.7 days with 
re-cahbrated model runs. 
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Figure 12: Parameter histograms after re-calibration with PEST posterior Monte Carlo 
analysis for test B8. Red line indicate PEST-estimated optimal parameter values and green 
lines are PEST-estimated 95% confidence intervals. 
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Figure 13: Posterior model parameter distributions estimated with DREAM for test B8. 
Red curves are normal distributions corresponding to mean and variance computed from 
data. 
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Figure 14: Model fits to test B8 data where only 30 data points were used in the optimization. 
The AICc for each model is included in parenthesis. 
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Figure 15: Temporal residuals of tests B4, B5 and B8 and the histogram of the test BJ 
residuals. 



26 



Table 1: 


PEST-Estimated Parameters Using Conservative 


Tracer Breakthrough Data 
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2.10 


610.1 


0.998 


13 






Single 


0.071 


14.3 








0.479 


0.989 


717 


B7 22Na 


0.5 


Double 


0.066 


9.89 


0.473 






0.601 


0.999 


183 






Multirate 
Single 


0.061 
0.068 


8.46 
12.9 


(0.398, 0.314) 


-0.921 


0.831 


0.580 
65.2 


0.999 
0.997 




496 


B8 22Na 


0.5 


Double 


0.065 


7.01 


0.289 






64.6 


0.999 


379 






Multirate 


0.047 


3.01 


(1.22, 45.8) 


-1.53 


1.86 


64.8 


1.000 









Single 


0.068 


15.8 








see a 


0.996 


676 


B{4,5,8} 




Double 

Multirate 


0.066 
0.045 


9.38 
3.46 


0.230 

(1.48, 107.9) 


-1.57 


1.98 


see h 
see c 


0.998 
0.998 


75 




a: tinj={305.8, 619, 
h: tij,j={306.6, 619, 
c: tinj={305.2, 615. 



3, 65.0} hours 
0, 65.3} hours 
2, 64.9} hours 



for B{4,5,8} 
for B{4,5,8} 
for B{4,5,8} 





Table 2 


PEST-Estimated Parameters Using 232ij Breakthrough Dat 


a 


Test 


Model 


-Rm 


^im 


ujd 


/i 


a 


tinj (hours) 


K' 




Single 


3.65 










1.92 


0.946 


B3 


Double 


2.36 


1.80 


0.754 






2.33 


0.995 




Multirate 


1.63 


5.68 


(1.44, 176.4) 


-1.86 


2.11 


3.00 


0.998 




Single 


3.49 










2.33 


0.987 


B7 


Double 


3.52 


2.83 


0.022 






2.15 


0.993 




Multirate 


3.48 


1.30 


(66.6, 2.45 X 10^) 


-2.21 


3.58 


2.12 


0.991 



27 



